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SUMMARY 


This report discusses the work partially supported under NASA 
Contract NAS 8-29662, 11 Segregation Effects During Solidification 
in Weightless Melts." The contract covers the period from July 5, 
1973 to July 4, 1974. 

During the contract period, the generalized problem of determin 
ing the temperature and solute concentration profiles during direc- 
tional solidification of binary alloys with surface evaporation has 
been mathematically formulated. Realistic initial and boundary 
conditions have been defined, and a computer program has been de- 
veloped and checked out. 

The program computes the positions of two moving (evaporation 
and solidification) boundaries and their velocities of movement, 
and also the temperature and solute concentration profiles in the 
semi-infinite material body at selected instants of time. 

The program has the following unique features: 

• Two moving boundaries are involved, i.e., 
the evaporative boundary and freezing 
boundary 

• Surface evaporation, and its related ef- 
fects such as material loss, evaporative 
segregation, and surface cooling due to 
the heat of evaporation, have been con- 
sidered 

• Surface temperature is realistically de- 
termined by the combined effect of heat 
radiation, evaporative cooling, and 
thermal diffusion 
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• Material parameters such as solid and liquid 
densities, specific heats, thermal conductivi- 
ties, mass diffusivities , and latent heat of 
fusion or evaporation, can all vary with both 
the temperature and composition 

• Realistic phase diagrams involving curved 
liquidus and solidus lines are used 

Our computer simulation work on solidification clearly shows 
that constitutional supercooling readily occurs and within-melt 
nucleation must then happen, particularly with reduced effective 
liquid mass transfer under zero-gravity conditions. Such results 
enabled us to explain and correlate some perplexing space solidi- 
fication phenomena observed on Skylab, e.g., E. McKannan*s weld 
(M551) and Prof. Adams’ braze (M552) results (see Monthly Progress 
Reports Nos. 10 and 11). Detailed and quantitative application of 
the results of this computer program, however, still awaits the 
gathering of pertinent crystal growth data. A final report is ex- 
pected to be written after these data are gathered and correlated. 
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INTRODUCTION 


Space processing is moving closer to reality. Bigger, better, 
and more uniform single crystals of important semiconductors and 
welds or brazes of improved properties have already been made in 
space, as reported in the Third Space Processing Symposium at 
Marshall Space Flight Center. Although processing of structural 
materials may certainly have a profit potential in the long range, 
it appears that the high cost per pound of single -crystal elec- 
tronic and optical materials makes these materials the most de- 
sirable contenders for immediate profitable returns from space 
processing. A selected single crystal study is, therefore, highly 
desirable to help us understand the segregation effects during 
solidification in weightless melts. 

Important tools for understanding these segregation effects 
are analytic solutions or computer programs that simulate or pre- 
dict what actually happens during space manufacturing. Such solu- 
tions and programs, furthermore, are probably necessary in space 

processing and other experiments where available time and experi- 

< 

mental facilities are limited, the cost per sample or experiment 
is very high, and yet only a limited total number of tests or test 
samples can be conducted. 

Theoretical predictions often greatly save time while compu- 
ter simulation saves cost. Specifically, analytic solutions and 
computer programs allow us to answer many questions during the 
planning or execution of space experiments on material solidifica- 
tion, such as learning 

• What phenomena are most important and what other 
phenomena are negligible 
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• Which influences are favorable to our understand- 
ing of weightless solidification and which are not 

• What conditions lead to optimal combination of 
the favorable influences or elimination of the 
unfavorable ones 

• What sample and processing conditions should be 
used 

• What is the best way to analyse the resultant 
samples for understanding a particular phenomenon 
or influence 

• How to save time and money — that is, how to 
maximize scientific return 

We have developed a number of analytic solutions relating to 
solidification and evaporation (Refs. 1-3). Several important com' 
puter programs have also been developed. Some of these solutions 
and programs were developed under our Contract NAS 8-27891, and 
they are already proving useful in correlating actual experimental 
results (Refs. 4 and 3). 

These analytic solutions and computer programs are, however, 
still in their early stages of development. The physical models 
involved are very simple and require considerable improvements to 
be used for other applications. It is, therefore, an important 
objective of this contract to refine and improve these models and 
the resultant analytic solutions and computer programs. 

These refined solutions and programs are more widely useful, 
have greater predictive value, and provide more accurate results. 
Such accuracies are absolutely necessary to separate the rather 
subtle zero-gravity effects on solidification, in the presence of 
noise due to other unavoidable or unanticipated but ever-present 
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miscellaneous effects. As a result of this continued work, more 
efficient space experiments and greater scientific returns appear 
possible. More meaningful solidification experiments and fuller 
utilization of the unique space environment may also result. 

The predicted results of our refined solutions and programs 
should, of course, first be checked with selected experiments. 
Another objective of this contract is, therefore, to design unique 
experiments to correlate the numerical results to actual solidifi- 
cation processes. This work is yet to be reported. 

Review of Previous Contract 

Under our NASA Contract NAS 8-27891, ’’Segregation Effects Dur- 
ing Solidification in Weightless Melts” (Ref. 3), two types of melt 
segregation effects were studied: evaporative segregation, or 

segregation due to surface evaporation, and freezing segregation, 
or segregation due to liquid-solid phase transformation. 

These segregation effects are closely related. In fact, evapo- 
rative segregation always precedes freezing segregation to some de- 
gree and must often be studied prior to performing meaningful solidi 
fication experiments. This is particularly true since evaporation 
may cause the melt composition, at least at the critical surface 
regions or layers, to be affected manyfold, often within seconds, 
so that at the surface region or layer the melting point and other 
thermophysical properties, nucleation characteristics, base for 
undercooling, and critical velocity to avoid constitutional super- 
cooling, may be completely unexpected. 

■ i 

To predict the segregation effects of solidification time and 
temperature and to correlate these predictions with actual experi- 
mental data, "normal evaporation equations” were developed (Refs. 1, 
4-6) . An evaporative congruent temperature (or equi-evaporative 
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temperature) was then defined and listed for various binary or 
ternary alloys . Knowing these congruent temperatures and the 
solute and solvent evaporating rates, one can predict the type 
(solute depletion or enrichment) and magnitude of compositional 
or constitutional changes on the critical melt surface. One ap- 
plication of this unique temperature is to explain, predict, or 
plan "anomalous" evaporative or constitutional melting (on cooling) 
or solidification (on heating) experiments. We then computed for 
a simple model the reactive jetting forces due to surface evapora- 
tion and, in particular, showed that these forces can be very sub- 
stantial on a differentially heated sample and may completely 
destroy the unique zero-gravity environment in space manufacturing 
(Ref. 7). In addition, these jetting forces may initiate surface 
deformation and vibration or other fluid disturbances, and may even 
produce some convection currents not normally anticipated. These 
studies also showed which sample materials are preferable, which 
should be avoided, and what impurities are harmful in producing ex- 
cessive jetting or effective as stabilizing influences. The rela- 
tionship between normal evaporation and normal freezing was then 
considered. Finally, applications of evaporation to space manufac- 
turing concerning material loss and dimensional control, composi- 
tional changes, evaporative purification, surface cooling, mate- 
rials standards, and freezing data interpretation were briefly de- 
scribed. 

In the area of segregation due to solidification, we explained 
in some detail the normal freezing process and its successful use 
in the semiconductor industry. Various constitutional diagrams 
demonstrated the desirability of using nonconstant segregation co- 
efficient techniques in metallurgical studies. We then stated the 
basic normal freezing differential equation, together with its 
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solutions for cases where the liquidus and solidus are quadratic, 
cubic, high-degree polynomial, and exponential functions of the 
melt temperature. The meaning of constant segregation coefficient 
was discussed, together with the associated errors due to curva- 
tures of the liquidus and solidus lines and the best value of 
constant segregation coefficient for a given solidification ex- 
periment. Numerical methods for computing the normal freezing 
behavior were then given. Finally, as an example, the steady 
state solidification of the Ni-Sn system under conditions of 
limited liquid diffusion and nonconstant segregation coefficients 
was described. This system was studied in the M553 experiment on 
Sky lab. 
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IMPORTANCE OF EVAPORATION 


Evaporation is important in space melting and solidification 
for the following reasons: 

• Significant evaporation of alloy components 
always occurs at high temperatures in space 
vacuum environments 

• High-temperature evaporation of alloys is gen- 
erally a neglected area of systematic research. 

Yet, unless the complete evaporative segregation 
behavior is understood and analyzed, solidifica- 
tion and its related segregation effects may not 
be properly studied because of ill-defined ini- 
tial conditions. Before the liquid alloy can be 
controllably solidified or even melted, there is 
invariably some surface evaporation to cause 
changes in composition, freezing temperature, 
supercooling characteristics, nucleation and 
growth morphology conditions, and the like 

• Controlled space evaporation probably most 
closely meets the requirements of our model of 
normal evaporation. We may thus be able to ob- 
tain material purity or evaporation standards, 

i 

thermal properties, or even such basic thermo- 
dynamic properties as heat of evaporation, ac- 
tivity coefficients, and sticking coefficients 
that are difficult or impossible to obtain on 
earth 
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• Evaporation is a much simpler process than freez- 
ing, since the former does not involve such com- 
plicated phenomena as nucleation, phase transforma- 
tion, and constitutional or nonconstitutional 
supercooling. Thus, in normal evaporation for 
specific geometries or alloy systems, we may 
ideally isolate and investigate such other phe- 
nomena as heat conduction or radiation, liquid 

or solid diffusion, fluid dynamics, and convec- 
tion currents. Exact knowledge of these phenomena 
is necessary to understand solidification 

• Evaporation causes surface cooling due to the 
heat of evaporation. This evaporative cooling 
effect is particularly important in low-melting 
materials (Ref. 8) 

• Different rates of evaporation at various sur- 
face regions give rise to unbalanced forces and 
momenta that may produce erratic or unwanted 
accelerations, surface distortions and vibra- 
tions, exceedingly large ’’equivalent gravities," 
and possibly new types of powerful convection 
currents in zero-gravity conditions 

• Evaporation may cause the surface composition 
of certain unwanted or unsuspected impurities 
to be increased a thousandfold or millionfold 
within seconds so that the layer's melting 
point and other thermo physical properties, 
nucleation characteristics, base for under- 
cooling, and critical velocity to avoid con- 
stitutional supercooling may be completely 
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unexpected. In fact, anomalous "constitutional" 
or evaporative melting on cooling, or solidifica- 
tion on heating, is possible because of surface 
evaporation. In addition, very large artificial 
gravities (e.g., 10 g) , strong fluid disturbances, 

or even new and significant convection currents may 
be produced from surface evaporation. These phenom- 
ena have been observed in the M553 movies, according 
to Dr. Martin Tobin of Westinghouse Co., Pa. 

The much greater evaporative segregation effects, if unac- 
counted for, would almost certainly conceal any minor or subtle 
zero-gravity effects, particularly in the presence of other unknown 
or uncontrolled effects. Definitive space solidification work 
should probably, therefore, be preceded by an evaporative compati- 
bility study of the sample materials and their possible associated 
impurities. In fact, evaporation is almost certain to be very im- 
portant or so overwhelming that the effect of zero-gravity or 
freezing segregation may be masked or even reversed. A freely 
suspended molten drop in space may, for example, have its surface 
solute concentration greatly enriched (as much as a millionfold), 
by neglected and undetectable trace impurities within seconds of 
its deployment. We are then dealing at the critical surface layer 
with a completely new and unanticipated alloy having an entirely 
different composition, melting point, surface tension, thermophysi- 
cal properties, latent heat of fusion, undercooling and nucleation 
characteristics, growth morphology, and the like. 

From this we can also see that any analytical, numerical, or 
experimental study on solidification may yield completely unex- 
pected or irrelevant results if the important and ever-present 
evaporation phenomena is not adequately taken care of. This is 
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particularly true in the study of nucleation, undercooling, and 
space manufacturing. Another important aspect of the present con- 
tractual work is to incorporate this generally neglected evapora- 
tion phenomena to define the exact initial and boundary conditions 
before and during the alloy solidification process. 
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COMPUTER PROGRAM WITHOUT SURFACE EVAPORATION AND RADIATION 


Solidification, even in one-g, is a complicated process in- 
volving a multitude of interrelated phenomena such as mass and 
heat transfer, phase change, and fluid motion. Comprehensive re- 
views on solidification have been given, for example, by Chalmers 
(Ref. 9), Tiller (Ref. 10), Christian (Ref. 11), and Li (Ref. 12). 

Solidification in zero-g is still very complicated. Here, 
gravitational force is negligibly small, but other effects as a 
result become important. For example, surface tension often 
plays a dominant role in determining the sample shape, processing 
technique and the resulting contamination level of the processed 
samples. Evaporation is another ever-present, complicating or 
dominating factor, but one that may be used to advantage when 
understood. Neglected, or improperly controlled evaporation may 
drastically change sample surface composition, fluid motion, equiva 
lent gravities, nucleation, and undercooling characteristics as 
previously described. The previous program, under Contract 
NAS 8-27891, however, does not deal with evaporation. 

Mathematical Definition of Solidification Proble m 

To understand thoroughly solute segregation either from com- 
bined evaporation and solidification, or in single-crystal growth, 
one requires a complete characterization of the (mass) diffusion 
and temperature fields in the solid crystal and remaining melt. 

The zero-gravity effect on the solidification may be overshadowed 
by other effects invariably present (such as evaporation) in any 
such growth process — a condition necessitating that such charac- 
terization be accurately defined. Unfortunately, the coupled par- 
tial differential equations for the diffusion and temperature 


13 



fields are generally not solvable. Although special case solutions 
have been given for some types of usually physically nonsatisfying, 
two-phase Stefan problems, for the general case solution we must 
resort to numerical computations. Existing numerical methods are 
always subject to such unrealistic assumptions as constancy of 
interfacial velocity, temperature or temperature gradients, segre- 
gation coefficients, diffusion constants, and other material thermo- 
physical properties. 

Under NAS 8-2,7891, a number of computer programs were de- 
veloped to study the unidirectional solidification of a binary 
alloy. These programs employ analytical and numerical methods. 

The analytic program is based on some closed-form solutions of a 
simple model and gives results for our numerical program to com- 
pare. The model for the analytic program deals with a binary 
alloy at a constant temperature and concentration throughout the 
initial liquid melt, with the surface temperature instantaneously 
dropped below the liquidus temperature. The liquid-solid inter- 
face temperature is assumed constant, and the concentrations of 
the alloy at the interface are given by the phase diagram having 
curved liquidus and solidus lines. In addition, the interface 
boundary plane moves according to a square root law relative to 
the solidification time. The program also allows the interface 
temperature and interface boundary to vary from these fixed rules, 
but in practice the variation is negligible and not above the com- 
puter error level (Ref. 3) . 

Although covered in detail in the final report on NAS 8-27891, 
the mathematical formulation of the model is presented below for 
the sake of completeness. 

We deal in unidirectional solidification with a liquid binary 
alloy to be directionally solidified into two phases, liquid and 
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solid. We consider the liquid alloy to be semi-infinite with origi 
nal (at t = 0) temperature T q and concentration C Q . Solidifi 
cation occurs when the temperature at x = 0 is changed from T q 
to a lower value T^, either instantaneously or gradually, so 
that T-^ is below the temperature T 2 at which the liquid mix- 
ture at concentration C q can be in equilibrium with a solid 
phase. As solidification occurs, the solid phase grows and its 
boundary is located at x = y(t) , and the interface temperature 
at this point is T\(t). The partial differential equations de- 
scribing the solidification process are the following: 


-,2 m 
0 T 

c 

St 

O. 

b 2 C 

Q 

bC 

0 

bx 

bt ’ 

D s > 2 

OX 

- St 



D ^ 


bx 2 

m 

bt * 

Sx 2 

_ 

bt 


for 


for 


0 < x < y(t) 


y (t) < x < 00 


( 1 ) 


( 2 ) 


where the variables T, C represent the temperature and concentra- 
tion (of solute in solvent) and the subscripts i ? s denote the 
liquid and solid phases, respectively. The thermal and mass diffu- 
sion coefficients a g , a D g , are assumed constant. The fol- 

lowing conditions are usually assumed throughout: 


(a) 

T^(x,0) = 

T and C,,(x,0) = C 

0 l ' 0 

(b) 

TgKt) = 

T and = C 

(c) 

T s(y(t),t^ 

1 = Tg(y(t) ,t) = T t (t) 

(d) 

c 8 (y ( t),t) 

' = f s ( T i (t) ) 
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(e) 

c^(y(t) ,t) 




• 

(f) 

PYy(t) = k s - ki for 

x = y (t) 

(g) 

- £ i( T i<t>)]y(t) = Dj 

n 

bx s 


for x = y(t) 

In many cases, it is also assumed 
00 y(t) = a . 

Equation (a) describes the condition that the original mix- 

ture is all liquid at temperature T and concentration C . 

° ° 

Equation (b) is a consequence of the semi-infinite nature of the 
mixture so that at any time t, the portion near infinity is un- 
changed. Equation (c) assumes that at the solid-liquid interface 
plane there is an interface temperature Th(t) and that both the 
solid and liquid phases at x = y(t) have this temperature. There 
is no discontinuity in temperature. Equations (d) and (e) state 
that the concentrations of solid and liquid at the interface are 
given by the solidus and liquidus curves, respectively, of the con- 
stitutional diagram for the alloy. Equation (f) connects the de- 
rivative of the moving boundary with the redistribution of tempera- 
ture and Eq. (g) connects the same boundary with that of concentra- 
tion. Equation (h) relates the position of the interface boundary 
to the solidification time t. 

The conditions on T (0,t) and C (0,t) are not fixed in 

o s 

our discussion, and a number of alternatives are considered: 
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1. T (0,t) = T-, (t) with T-, (t) equal to a constant for 

S -L 1 

all t ; 

2. linear, T., (t) = T + t(T-, - T )/s for t < s and 

I ' o v 1 o 

T^t) = T l for t > s; 

3. exponential, T^(t) - + (T - T^) e t ^ S so 

T l(0) = T q and T x (») = t 

For C (0,t) the conditions considered are C (0,t) = C 1 
s si. 

usually taken CgC^) or at times a condition conserving mass 
between 0 and °o. 

The two approaches we have pursued may be designated as ana- 
lytic and numerical. The numerical approach can be applied to all 
three conditions on temperature whereas the analytic approach holds 
only the case of constant temperature instantaneously applied. A 
variant of this analytic method to apply to linear varying tem- 
perature has been investigated. 

An analytic solution to the coupled partial differential equa- 
tions (1) and (2) subject to the initial and boundary conditions 
(a) through (g) has been given (Ref. 13). A numerical program has 
been designed for the analytic solution. 

These numerical programs developed under NAS 8-27891 are 
based upon finite difference approximations of the partial and 
ordinary derivatives and involve a variable spacing (for improved 
computing efficiency) . The programs have given acceptable results 
and compared well with the reference analytic solution, where com- 
parable. The basic physical properties such as densities, diffu- 
sivities, specific heats, thermal conductivities, and heat of 
fusion have been held to be constant, and independent of tempera- 
tures and concentrations. 
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COMPUTER PROGRAM DEVELOPED UNDER PRESENT CONTRACT 

Under the present contract, we have extended the programs to 
allow for reasonable variation of these physical properties. The 
approach that has first been taken is to base the values of these 
physical properties upon extrapolated values of temperature and 
concentration, and then to determine the values of temperature and 
concentration. The process is then repeated by re-evaluations of 
the physical properties. Other modifications of our original pro- 
gram are: 1) to store the physical properties for each of the 

mesh points and to employ the appropriate quantities at each step, 
and 2) to recheck the mass and heat diffusion equations to make 
certain that the constancy of these properties is not assumed. 

An additional major program modification has been the inclu- 
sion of evaporation effects. This includes evaporation before 
solidification that is mathematically identical to the problem of 
simple solidification in binary alloys. After solidification 
starts, significant evaporation may still exist. We then have to 
deal with two moving (solid-gas and solid-liquid) boundaries 
located at y(t) for evaporation and at z (t) for solidification, 
as will be described. 

Modification of the initial and boundary conditions a-h has 

also been made to make the problem more physically meaningful. One 

such modification is to include a surface heat radiative loss term 
4 • 

involving T . This term affects the convergence of the problem 
and creates the need for different algorithms. As reported pre- 
viously, (Ref. 14), the surface cooling due to evaporation is neg- 
ligible for many metallic systems such as nickel and iron alloys, 
or other higher melting materials, and has not been studied in de- 
tail under this contract. 

wuemma fid bum m man 
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To obtain solutions for realistic boundary conditions and to 
include the various mass transfer effects, numerical solutions of 
the partial differential equations of heat and mass transfer are 
required. We have again used the finite difference method to ob- 
tain the numerical solution. 


The boundary conditions for surface temperature include radia 
tion cooling as given by the Stefan-Boltzmann equation and also 
include evaporative cooling for both components of the alloy. 

D JlAi 1 1 f" 1 O 1 Qt.t 1^ qq noci jA -i -rs vs -? 4-1* ~ J 

nu,u wwtu u.4juuiiijcu. UC LCllUJ.UjLi.lg LUC ti V <X]JKJ L Ct U J-Uli Ld-LCis 

At the interface it is assumed that the temperature and concentra- 
tion relationships for each phase are given by the constitutional 
diagram for the alloy. The temperature dependence of the thermal 
and mass diffusion coefficients are allowed for each phase. 
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GENERALIZED SOLIDIFICATION WITH SURFACE EVAPORATION 


Evaporative Solidification of a Binary Alloy 

Given a semi-inf inite binary alloy melt, initially at concen- 
tration C^ and temperature T q , we consider the solidification 
of the alloy due to surface heat loss by evaporation and radiation 
(Fig. 1). There are two separate regimes to be considered. The 
first is concerned with temperature and concentration variations 
before solidification begins; the evaporation causes the original 
liquid-vapor boundary to change. Thus, we have a moving boundary 
problem. The second regime begins with the solidification which 
introduces a boundary between the freezing solid and remaining 
liquid phases whose compositions, we assume, follow the phase 
diagram, i.e., solidus and liquidus curve relations hold. Conse- 
quently, after solidification begins, there are two moving bounda- 
ries: one is the evaporative boundary and the other is the freez- 

ing or solidification boundary. 



Fig. 1 Evapor at ion -Solidification of 0.01 Sb-Ge 
Initially at 970 °C at 0.01 Second 
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Equations at the Evaporative Bounda ry 


We denote the evaporative boundary as x = y(t) where 

2 

y(0) = 0. The evaporation rates in mol/m /sec for pure solute 
and solvent are, respectively (Ref. 15), 
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A -B / T 
u u 


(M T ) 
u s 7 


1 

2 


A _T> /»T 

in. u f L 

V = K 10 v v (M T )" 2 
e vs 7 

-5 

where K = 5.83 x 10 , M , M are molecular weights for solute 

G II V 

and solvent atoms, T is the evaporating surface temperature in 
degree K, and A^, B^, A^ t B v are the evaporating constants for 
solute and solvent, respectively. If and p^ are the solute 

and solvent densities, then 

dy = VM v (l - C) 

dt p + p 

r u r v 

where C is the concentration at the moving boundary. 

The heat loss rate equation at the boundary due to radiation 
and evaporation is given by 

If = - €oT s - U V - - c > 

where e is emissivity coefficient, a the Stefan-Boltzmann con- 
stant, and 7^ and y^ are specific heats for solute and solvent, 
respectively. 

The equation for the rate of concentration change is 



Since the evaporative boundary is a moving one, and since both the 
evaporation temperature T and solute concentration C are func- 
tions of distance x *= y(t) and time t, i.e., T = T(x,t) and 
C - C(x,t), the total derivatives and may be obtained 

from the partials, i.e.. 


dT _ ST /dT\ 

dt St \Sx/ \dt/ 

x== y 


dc = sc + (h£\ (*2S 

dt St \Sx/ \dt/ 

x=y 


St d c 

where and ^ are evaluated at the moving boundaries . 

dy ht dC 

Given — , -jp and — , we can integrate for y, T, and C 
for the moving boundary using a modified Euler method. 


v^ = v + At 
t+At 


(i+1) . At 

V t+At = V t + ~ 


/dv\ 

/"dv' ) 

VdJ + 

\dt/ 


(i) 

t+At 


where 


dv 


(i) 


dt 


using the value v 


is the value of the derivative at time t + At 

(i) 


for v. 


ST 


To determine — at time t and t + At requires knowledge 
of the distribution of temperatures at both times. Those at time 
t+At are initially approximated by an extrapolation and are cor- 
rected using an approximated value of the temperature of the evap- 
orating boundary with the heat diffusion difference equations. 
Since the change in temperature at the boundary is greatest due to 
the heat of evaporation, more iterations are applied to determine 
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it than to the determination of temperature distribution by means 
of diffusion equations. Similar considerations hold for the de- 
termination of g— and 

The computations of the position of the evaporation boundary 
ty - y(t) ], temperature (T) , and solute concentration (c) at 
this boundary constitute an initial value problem in ordinary dif- 
ferential equations. Thus, given y =0, T=T,C=C, at 

o o o 

time t = 0, and given also the equations for velocity of movement 
of this boundary dy/dt, and rate of change of temperature and 
solute concentration dT/dt, and dC/dt, we can determine for 
selected times the values of y, T, and C. The method used is 
an iterated Euler scheme: 

y n+l = y n + f (y n + y n+l } 

where the initial value y' is taken as y'. This scheme must 
be connected to the problem of determining the temperature and 
solute concentration distribution within the semi-infinite body 
because the derivatives dy/dt, DT/dt, and dC/dt depend upon 
these quantities. The first step is to determine a first approxi- 
mation of the temperature and solute concentration by extrapola- 
tion and then correct these values from the newly approximated 
values of the boundary position and the temperature and concentra- 
tion thereat. 

Start of Solidification 

To determine the time when solidification has begun, the 
boundary temperature is compared with the temperature obtained by 
the inverse function for the liquidus curve evaluated at the 
boundary concentration. If the former is greater, then solidifica- 
tion has not yet begun. If it is smaller, then solidification has 
begun. In order to avoid an exact iterative procedure to determine 
the instant of solidification and to follow it up by a starting 
procedure for the first time interval thereafter, a simplified 
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approach has been taken that introduces a small error in the evapo- 
rative boundary and freezing boundary. By allowing the temperature 
to be below the solidification temperature by a small amount and 
by assuming that the temperatures at both boundaries are the same, 
a starting value of x = z(t) of the freezing boundary is deter- 
mined so that the loss in concentration due to solidification is 
compensated by the gain in concentration at the liquidus. Given 
the new temperature TI2 below the temperature at which solidifi- 
cation begins, we compute CSS = FS(TI2) and CLL = FL(TI2) , the 
corresponding solid and liquid concentrations given by the phase 
diagram. To determine DEL2 - ZI2 - YI2, the distance between 
the evaporative boundary and solid-liquid interface, we assume 
that the solid is entirely at concentration CSS, and the liquid 
varies linearly from CLL to CC(II2), the concentration at the 
first mesh point x(II2) after the evaporative boundary. The 
total concentration is to equal the concentration in the whole 
regime had no solidification taken place. We assume it to be 
CL2 computed at YI2 and to vary linearly to CC(II2) at 
x(II 2). This yields the equation 

CSS * DELZ + (cLL + CC(II2))/2 * (x(II2) - YI2 - DELz) 

= (cL2 + CC(H2))/2 * (x (112) - YI2) 

Hence 

DELZ * [CSS - ( CLL + CC(II2) /2) ] = (x(II2) - YI2) * (CL2-CLL)/2 
where 

DELZ = (CLL - CL2) /2 * (x(II2) - YI2) -s- (cLL+ CC(II2)/2 - CSs) 


Then 


ZI2 = YI2 + DELZ 


and 


dz _ DELZ 
dt BELTS * 
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This enables us to begin the next time step with initial 
values for y(t ), z(t ), jjp #, l(y(t )) - T(z(t )) = TI2, and 

c(y(t s )) = css, c s ( z(t ;)) = css! Cl (.(t;>).- CLL. 


The Two -Boundary Problem" Derivative Estimation 


The equations at the freezing boundary are those given in 
the Grumman Final Report RE-458 to Contract NAS 8-27891 (Ref. 3) , 
with the exception that the freezing boundary is now called 
x = z (t) and not x = y(t) as in Eq« 49 c-^g. At every time step 
we must compute (in addition to the temperature and concentration 
at the evaporation boundary) the temperature at the freezing bound- 
ary. The concentrations are determined by the phase diagram. The 
method we employ is that which determines T (the solidification 
temperature) and by means of Eq. (49) f,g, Having obtained 

we obtain z(t) by means of a modified Euler method. Since 
the Eq. (49) f,g required approximation for (g^ S ’) 2 t and (§x^) t* 

we must develop techniques for these approximations appropriate to 

various situations for mesh points. In addition, for the computa- 
. dy dT dC 

tion of and — at the evaporative boundary, we also 

need (tx)y, t and (§x)y,t* When there are two mesh points be- 
tween y and z, then the techniques alluded to above are avail- 


.2 

able. This involves determining — ~ at both y and z 

9 dx 

o C 

same for «*. When there is only one mesh point between 
dx 
S^T 

z, then t: at both points are the same. When there are 

Sx Z 

points between y and z, then we can assume either that 
is zero and hence (|g = (g) z = ~ or that 


and the 
y and 
no mesh 
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s 2 t 

9x 2 


= k and hence (^~) “ (l - 


z - y 


■)S 


(IS, ■ (i - 


z - y 


\ (3T 

/§ 3 C* The c ^°^ ce °f k xmst be small so that 


s 2 t 

Sx 2 


= k 


5 2 T 

Sx 2 


= k 


2 St 
S x 


is negligible. Thus, since z-y is also 

v 2 

S T 


very small this option is indistinguishable from — y =0. We 

Sx 

have three cases: 1) no mesh points between two boundaries and we 

S^t St S 

assume — j ~ (sx) ” (sSz* one mes ^ point between y and 
Sx x ^ z 

S^T 

z when — is obtained from the three points and 
Sx Z 


/ST\ _ T (z) - T(y) (z - y) S 2 T 


\Sx/y 


z-y 


and 


/ST\ ^ TOO ~ T(y) 


Sx' 


\Sx/ 


Z-y 


(z - y) S T 

Y~ 2 *, and 3) when two or more mesh points, say x. and 


Sx' 


x 


/b^T\ 

are between y and z so that we can compute ^ and 

/ ST\ _ T(Xj) - T(y) 


/S 2 T\ 

\ — 2/ separately and distinct. Then 

Sx y 


VSx/ y 


x. - y 

i J 


(x. - y) 

1. | 

6 

CM 

ro 

i and I 

y 

f dT\ 

T(z) - T(x. +1 ) ( (z - x. +1 ) 

fS T\ 

2 1 

U 2) 

Sx 


z 2 - x i+l 2 

Sx^ 


S^T S^C 

In general, it is necessary to compute — y and — in three 

Sx Sx 

ways, two ways indicated above for the solid regime and a third 
for the liquid side of the freezing boundary. It is similarly 
necessary to compute and in three ways. 


B oundary and Mesh Points 

When boundary points come close to mesh points, the computa 
tion of derivatives may be vitiated by closeness to mesh point. 
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Therefore, tests are made to determine when such closeness occurs 
as usually expressed in terms of a decimal fraction of the interval. 
In that case, the reference point is moved to the next mesh point 
and the values of T and C at the skipped mesh point are ob- 
tained by linear interpolation. This interpolation depends on 
which side of the solid-liquid interface the mesh point lies. For 
the evaporative boundary similar considerations hold. 

Solution for R emaining Poin ts 

The solution for the remaining points is obtained as in the 
Final Report previously mentioned, pages 3-14 and 3-15 (Ref. 3). One 
change is, however, necessary because the first mesh point (or 
more) are no longer under consideration if the evaporative boundary 
has passed them. The subroutine TRIST is used to solve for the re- 
maining points. In this subroutine we compute the values of tem- 
perature and concentration at intermediate mesh points when given 
the values at the two extreme mesh points. We replace the values 
at the mesh point to the left of the evaporative boundary by those 
at the evaporative boundary point, before solving for the inter- 
mediate points. This can be done without destroying any useful 
information since that mesh point is no longer used in the compu- 
tations. The subroutine TRIST does not depend upon equal spacing 
or any regular spacing and therefore can accommodate this usage. 

Conv ergenc e 

The convergence problem is the crux of the program. Oscilla- 
tion tends to cause the needed quantities to overflow. Thus, tests 
must be made on all the quantities to contain them within reason- 
able bounds. The subroutine MOTON is used to check the monotonicity 
of these consecutive points. In addition, the temperature at the 
evaporative boundary is necessarily less than the temperature at 
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the freezing boundary. This condition is always imposed in the 
program. 

In addition, the solution for the solidification temperature 
and freezing boundary derivative (especially the latter) involves 
very rapidly changing quantities. More iterations should, there- 
fore, be expended in this part of the program. Fewer iterations 
are needed for determining the evaporative boundary, and the tem- 
perature and concentration at that boundary. The program allows 
five iterations in the former for each of the latter. The number 
of iterations of the latter is used in a manner analogous to that 
described in Final Report RE-458 (Ref. 3) . 

An input quantity NIT (usually a multiple of 4) gives the 
maximum number of iterations. When NIT/2 iterations occur and 
convergence is not reached, the time step size is halved. This 
process is continued until either convergence is attained or the 

minimum step allowed by the program has been iterated NIT 4- 1 

0 

times. In this case the program may stop or continue on using the 
nonconverged quantities. Very often these quantities are suffi- 
ciently smooth so that convergence will occur on the next interval 
and the program gives satisfactory results. 

However, if the program proceeds with the minimum step and 
the maximum number of iterations, the results may be spurious. In 
case of overflow, there is no doubt of it. Otherwise the user 
must look at results to decide whether he finds them reasonable. 


c 


( 


29 



IMPROVED COMPUTER PROGRAM 


The complete computer program for the generalized solidifica- 
tion problem is listed herein (see appendix) , together with a 
glossary explaining the special names used in the program. This 
computer program has the following unique features: 

• Surface evaporation, and its related effects 
such as material loss, evaporative segrega- 
tion, and surface cooling due to the heat of 
evaporation, have been considered 

• Material parameters such as solid and liquid 
densities, specific heats, thermal conductivi- 
ties, mass diffusivities , and latent heat of 
fusion or evaporation, are allowed to vary 
with the temperature and composition 

• Realistic phase diagrams involving curved 
liquids and solidus lines are used 

• Two moving boundaries are involved, i.e,, 
the evaporative boundary and freezing 
boundary 

• Surface temperature is determined by the 
combined effect of heat radiation, evapora- 
tive cooling, and thermal diffusion 

Use of Computer Program 

The computer program works well if the following three input 
program parameters are properly chosen: 1) time step size (DELT) , 

2) grid spacing (DELX) , and 3) maximum iteration count (NIT) . 




^ nor 
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The solidification boundary is sensitive to the grid spacing. 
This is because in passing through a mesh point, discontinuity in 
the computation occurs for the following reasons. We compute 
the derivatives in terms of the temperatures and concentrations 
at the discrete mesh points. When one mesh point is dropped be- 
cause solidification occurs near it, the derivative based on a 
substituting new mesh point is discontinuous with that based on 
the previous mesh point. Though this discontinuity can be reduced 
by using a smaller time step size, it would be a self-defeating 
strategy. An alternative is to accept the discontinuous results 
as they occur, advantages being taken of the fact that the program 
corrects itself. Although the derivative dz/dt is large when 
the solid-liquid interface passes through a mesh point, it becomes 
smaller thereafter thereby correcting the solidification boundary 
position. 

The frequency of this self -correction depends on the grid 
spacing. Too small a grid spacing would cause too frequent self- 
corrections. Too large a grid spacing, on the other hand, would 
obscure the rapid temperature variations around the solidification 
boundary. This indicates that a proper choice of the grid spacing 
is required to achieve an optimal tradeoff between accuracy and 
computing time. There is another tradeoff between the time step 
size and maximum iteration count for optimal computing results. 

Since each evaporation-solidification problem represents a 
different and unique physical situation, each case must be dealt 
with separately. However, based on our experience, the following 
guidelines would be helpful: 

The first consideration for the choice of the grid spacing 
is the behavior of the evaporation boundary after solidification 
begins. If the evaporation boundary is virtually stationary as 
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compared to the solidification boundary, the grid spacing should 
be chosen so that the evaporation boundary is within the first 
mesh interval (between the first and second mesh points) . If, on 
the other hand, the evaporation boundary is moving at velocities 
comparable to those of the solidification boundary, then the grid 
spacing can be selected more freely. The major consideration in 
this case is the relationship between the grid spacing and the 
time step size. For a fixed time step size, the grid spacing 
should be chosen so that at least four time intervals (of step 
sizes) occur before the solidification boundary passes through a 
mesh point. 

In cases where the evaporation boundary is virtually sta- 
tionary, one must experiment to determine an optimum time step 
size in terms of accuracy and computing time. The conditions of 
the experiment are as follows. Set both the minimum time step 
size (DELTM) and the time printing interval (DELP) to zero. Set- 
ting the time printing interval to zero will cause the computer to 
print out every computer time step. Setting the minimum time step 
size to zero will not cause the program to cut back indefinitely 
but will use, as the minimum, the time step size divided by 1024. 
By examining the computed results, one can see at what time step 
sizes the program is running. By examining the actual iteration 
count (IT), one can see if the program is converging or not. If 
not converging repeatedly, a smaller time step size is indicated. 
If the program is converging most of the time, then the minimum 
time step size can be set at the level of the most frequent time 
step size and the actual iteration count re-examined to see if the 
program still converges most of the time. For long runs, the time 
printing interval, must not be zero or small, but must be chosen in 
consideration with the amount of the required output. 
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To improve the computing time on long runs, one should con- 
sider enlarging the grid spacing as suggested above as one of the 
tradeoffs. In addition, one may change the maximum iteration 
count upwards or downwards to also improve the computing time. 

Our computer program has the capability for assuming equal or 
unequal (doubling)- mesh point spacings. Our experience, as indi- 
cated in Tables 1-4, shows that the unequal spacing scheme gives 
practically the same accuracy with far less computations as com- 
pared with the equal spacing scheme. This may be due to the 
rapidity at which the temperature declines at the evaporation 
boundary. Other physical situations may give different results and 
may indicate that the equal spacing scheme should be used. 

The program input parameters consist of a set of integers IX, 
IAM, NIT, IM, N, and NCN; and a set of real numbers DELX, DELT, 
DELTM, DELP , TF , and S. IX is the maximum number of mesh points 
to be used in the program. Present, IX 28. IAM is the spacing 
option indicator. If IAM equals 0, the points of mesh are equally 
spaced with grid spacing DELX. If IAM = 1, an unequal spacing is 
indicated. The first two intervals are equal and set to DELX. 
Thereafter, each interval is double the previous interval in spac- 
ing. NIT is the maximum number of iterations as interpreted in 
the context of halving the time step size. If the step is begun 
at the minimum time step, the NIT is the maximum number of itera- 
tions allowed. IM is the number of mesh points in actual use. 

Ihe input value of IM introduces the minimum number of mesh points 
to be used. Thereafter additional mesh points are added as re- 
quired by a substantial change in temperature at next to last mesh 
point, that is, 1 degree below the initial temperature. IM is 
increased until IM is equal to IX. 
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TABLE 1 VARIATION OF TEMPERATURE (°C) AT EVAPORATIVE BOUNDARY 

FOR VARIOUS COMPUTATION SCHEMES 


Scheme 

I 

II 

III 

IV 

V 

Grid 

0.01 

0.01 

0.001 

0.001 

0.0001 cr 

Spacing 

equal 

unequal 

equal 

unequal 

unequal 

time, ms 

0.2 

966.5 

966.5 

966.5 

966.5 

966.5 

0.6 

959.4 

959.4 

959.4 

959.4 

959.5 

1.4 

945.7 

945.7 

945.7 

945.7 

945.8 

1.8 

938.9 

938.9 

938.9 

938.9 

939.0 

2.0 

935.5 

935.5 

935.5 

935.5 

935.6 

2.05 

934.7 

934.7 

934.7 

934.7 

N.C. 

2.075 

934.2 

934.2 

934.2 

934.2 

N.C. 

2.0875 

934.0 

934.0 

934.0 

934.0 

N.C. 

2.09375 

933.9 

933.9 

933.9 

933 . 9 

N.C. 

2.1 

933,8 

933.8 

933.8 

933.8 

933.9 

2.1125 

933.6 

933.6 

933.6 

933.6 

933.7* 

2.1375 

933.2 

933.2 

933.2 

933.2 

933.3* 

2.1875 

932.4 

932.4 

932.4 

932.4* 

932.5* 

2.2875 

930.8 

930.8 

930.8 

930 . 8* 

930.9* 

2.4875 

927.7 

f 927.7 

927.2 

927.7* 

927.7* 

2.8875 

921.5 

921.5 

921.5 

921.5* 

921.5* 

3.6875 

5.2875 

909.3 

885.8 

909.3 

885.8 

909.3 

885.8 

909.3* 

909.3* 


“ft 

Hand interpolations 
N.C. not computed 
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TABLE 2 VARIATION OF POSITION (pm) OF EVAPORATIVE BOUNDARY 
WITH TIME FOR VARIOUS COMPUTATION SCHEMES 


Scheme 

I 

II 

III 

IV 

V 

Grid 

0.01 

0.01 

0.001 

0.001 

0.0001 cm 

Spacing 

equal 

unequal 

equal 

unequal 

unequal 

t ime 5 ms 

0.2 

0.122 

0.122 

0.122 

0.122 

0.122 

0.6 

0.351 

0.351 

0.351 

0.351 

0.351 

1.4 

0.752 

0.752 

0.752 

0.752 

0.752 

1.8 

0.927 

0.927 

0.927 

0.927 

0.928 

2.0 

1.009 

1.009 

1.009 

1.009 

1.010 

2.05 

1.029 

1.029 

1.029 

1.029 

N.C. 

2.075 

1.039 

1.039 

1.039 

1.039 

N.C. 

2.0875 

1.044 

1.044 

1.044 

1.044 

N.C. 

2.09375 

1.046 

1.046 

1.047 

1.047 

N.C. 

2.1 

1.047 

1.047 

1.047 

1.047 

1.050 

2.1875 

1.048 

1.048 

1.048 

1.048* 

1.052 

2.2875 

1.049 

1.049 

1.049 

1.048* 

1.053 

2.4875 

1.051 

1.051 

1.051 

1.050* 

1.055 

2.8875 

1.055 

1.055 

1.055 

1.054* 

1.059 

3.6875 

1.061 

1.061 

1.062 

1.061 

1.066 

5.2875 

1.072 

1.072 

1.072 

_ 

.. 


Hand interpolations 
N.C. not computed 
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TABLE 3 

VARIATION OF 

POSITION 

(pirn) OF SOLID- 

-LIQUID 

INTERFACE 

Scheme 

I 

II 

III 

IV 

V 

Grid 

0.01 

0.01 

0.001 

0.001 

0.0001 cm 

Spacing 

equal 

unequal 

equal 

unequal 

unequal 

time, ms 

0.21 

0.205 

0.204 

0.109 

0.109 

0.106 

0.24875 

0.211 

0.211 

0.179 

0.180 

0.408 

0.28875 

0.229 

0.229 

0.350 

0.350* 

1.03 

0.36875 

0.283 

0.283 

0 . 866 

0.860* 

2.76 

0.52875 

0.429 

0.429 

2.291 

N.C. 

N.C. 


TABLE 4 

VARIATION OF 

TEMPERATURE 

(°C) AT SOLID 

-LIQUID 

INTERFACE 

Scheme 

I 

II 

III 

IV 

V 

Grid 

0,01 

0.01 

0.001 

0.001 

0.0001 cm 

Spacing 

equal 

unequal 

equal 

unequal 

unequal 

time, ms 

0.21 . 

933.8 

933.8 

933.8 

933 . 8 

933.9 

0.24875 

927.7 

927.7 

927.7 

927.7* 

927.9* 

0.28875 

924.5 

921.4 

921.5 

921.4* 

921.6* 

0.36875 

0.52875 

909.3 

885.8 

909.3 

885.8 

909.3 

885.8 

909.3* 

909.4* 


Hand interpolations 
N.C. not computed 
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NONCN is a nonconvergence option. Failure to converge occa- 
sionally is not necessarily an indication of unacceptable results. 
Therefore, it is desirable to continue computations and examine 
the results to see if they are acceptable. This is done by setting 
NONCN to 1. If NONCN is set at 0, the nonconvergent results are 
not printed unless called for by the print interval. If NONCN is 
-1, the program stops on nonconvergent results. 

The quantity DELX is the grid spacing. Equal spacing and 
unequal double spacing both make use of this quantity as indicated 
in the discussion of I AM. The quantity DELT is the maximum time 
interval (step size) for computation. The quantity DELTM is the 
input minimum time interval. The program uses as its actual mini- 
mum the larger of the quantities DELT/ 1024 and DELTM. Thus, even 
if DELTM is set at 0, the number of halving on cutting back the 
time interval is at most 10 (2 10 = 1024) . The program in its pre- 
solidification phase starts with its actual time step DELTS set to 
DELT/8 and allows it to build up to DELT by quick convergence. 

On the other hand, near the beginning of solidification, 

DELTS is allowed to cut back to as small as DELT/ 2 56 in order to 
find an acceptable start of solidification. After solidification 
has begun, then the restriction of DELTS is between DELTK and DELT. 
If halving reduces DELTS below DELTK, it is set to DELTK. The 
quantity DELP is the present interval. If DELP = 0, then every 
time step is printed. TPR is the time for outputing results. TPR 
is set originally to DELP and after printout is reset to TPR + DELP. 
The program prints results if the time TIME1 at the end of the 
time step equals or exceeds TPR. The program does not attempt to 
set DELTS so that TPR — TIMEl. This is only a slight inconve- 

nience when the print interval is large as compared to DELTM. Gen- 
erally, we would like DELTM to be set close to the most frequent 
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DELTS provided that failure to converge does not ensue on a regular 
basis. TF is the final time of program. This means that if TIME1 
equals or exceeds TF, no additional time steps are taken. 

The decimal quantity S between 0 and 0.5 is used to de- 
termine closeness to a mesh point. If the boundary point (either 
evaporation or solidification) is such that it exceeds the point 
that divided the mesh interval surrounding the boundary point in 
the ratio (1-S)/S, then the mesh reference point for computation 
is moved to the next mesh point. The introduction of S is to make 
the transition due to passing a mesh point less abruptly discontinu- 
ous. The best values of S are between 0.05 and 0.15. For com- 
putations on the solid side of the solidification boundary, we con- 
tinue to use the old mesh points until the boundary point passes 
the point that divides the new mesh interval about the solidifica- 
tion point in the ratio S/(1-^S). This strategy causes a gradual 
transition from one mesh point to another. The integers III, 112 
are used as reference point indicators for the solid and liquid 
sides, respectively. For the evaporation boundary, 113 is used to 
indicate which points are used. 114 is used only to indicate the 
first mesh point to the left of the evaporation boundary. 


Typical Computer Input 

The definitions of the various inputs fed into the computer 
are given in the Glossary of Program Parameters. Typical input 
values are as follows: 


IX 

IAM 

NIT 

IM 

NONCN 


28 - maximum number of mesh points 

1, unequal, doubled grid spacing 

20, maximum number of iterations 

16, actual number of points in mesh 

0, allowing the program to continue when non- 
convergence occurs with no special printout 
of these results. 
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The alloy phase diagrams are determined from the five con- 
stants ET, EA, EB, EC, and ED, which define the liquidus C 0 and 
solidus lines C as two functions of the temperature, T: 

O 

C»(T) = ED x (ET - T) 2 + EC x (ET - T) 

C_ (T) = EB x (ET - T) 2 + EA x (ET - T) 
s 

In our example of 10 mole percent (C q = 0.10) of antimony 
in germanium initially uniform at 970 °C (T q = 970) 

ET = melting point of pure germanium = 956 °C 
EA = 0.128812 x 10" 3 

EB = -0.82218 x 10 -7 

EC = 0.466678 x 10 -2 

ED = -0.60466 x 10 " 5 

The evaporation constants for the solvent and solute as de- 
fined previously under "The Equation at the Evaporative Boundary" 
are: 

AU = A = 0.1115 x 10 2 
BU = B = 0.863 x 10 4 

u 3 

EMU = M = 0.2435 x 10 

AV = A = 0.1171 x 10 2 

BV = B = 0.1803 x 10 5 

EMV = = 0.726 x 10 2 

EK = K = 5.833 x 10 -4 
e 

The diffusion coefficient of the solute in solvent in the 
solid and liquid states are, respectively 

DS = D g - 0.10 x 10 -6 

DL = D| = 0.10 x 10" 3 
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The density p, and latent heat of evaporation, 7, of the 
pure solvent are, respectively, 

RHO = p - 5.32 

v 

GAMMA = Y v - 160 

Corresponding values for pure solute are: 

RHOU - p =6.68 

u 

GAMMAU = 7 “ 39 

u 

The above give two derived quantities: 

AL3 = ^ = k 4 /p v c 

ASS = a = k /p c 

S S V 

where 

CEE = c = 0.740 x 10 ^ = specific heat 

The two input parameters in the surface radiation terms are: 
EE = e .= 0.55 = emissivity coefficient, and 

-7 

SIG = a = 0.136 x 10 - Stefan-Boltzmann constant. 


Computer Output 

The first line of computer outputs gives the program input 
parameters IX, IAM, NIT, IM, and NONCN, which are defined pre- 
viously and also in the "Glossary." The next two lines of computer 
output give the phase diagram constants (ET, EA, EB, EC, and ED) 
and evaporation constants (AU, BU, EMU, A V, BV, EMV, and EK) , re- 
spectively. The next printouts are for CEE, DS, DL, TO, CO, XKL, 
RHO, GAMMA, RHOU, GAMMAU, EE, SIG, T2, and COO, where T2 is the 
temperature when solidification begins for the melt of initial 
solute concentration CO. 
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The computed numbers are then outputed as follows: 

IT 38 actual iteration count 

IM ■* number of points in mesh 

111 * grid point reference for solid side of mesh 

112 = grid point reference for liquid side of mesh 

113 = grid point reference for evaporation boundary 

114 = grid point reference for point after evaporation 

boundary 

These printouts are then followed by the computed values asso- 
ciated with the evaporation boundary: time, location y, concen- 

tration C, temperature T, extent of points X(IM) , current 
time interval DELTS, dy/dt DYDTl, dc/dt DCDTl, dT/dt DTDTl. If 
solidification has not begun, then there is no information about 
the solidification boundary. Otherwise, we have position of the 
solidification boundary z computer language (ZI1), solid solute 
concentration C g (CS1) , liquid solute concentration C^(CLl), 
temperature T (Til), and rate of movement dz/dt (DZDT1) . All 
decimal outputs are five per line with excess going to the next 
lines . 


Representative Computed Results 


The study of the effect of varying the grid spacing DELX on 
the computed results is summarized in Tables 1 through 4. The 
five cases considered are: 


Case I: DELX =* 0.01 cm with equal spacing 

Case II: DELX = 0.01 cm with unequal spacing 

Case III: DELX = 0.001 cm with equal spacing 

Case IV: DELX - 0.001 cm with unequal spacing 

Case V: DELX - 0.0001 cm with unequal spacing 
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Tables 1 and 2 indicate the insensitivity of the evaporation 
boundary and its temperature to grid spacing, provided that the 
spacing is always larger than the evaporation boundary point. 

Tables 3 and 4 involve solidification boundary and show that the 
temperature is insensitive to DELX but that the solidification 
boundary is quite sensitive to the choice of DELX. Thus, it is 
important to use DELX sufficiently small so that the solidifica- 
tion boundary movement is fully exhibited and not stunted by a 
large grid spacing relative to which the boundary size is small. 

The spacing affects the evaluation of the first and second temper- 
ature partial derivatives with time which are larger in absolute 
values for smaller spacings, due to more rapid temperature changes 
near the boundaries . 

The figures (Figs. 2-4), prepared from the computed results in 
Tables 1-4, indicate the superiority of unequal over equal grid 
spacing. For DELX = 0.01 cm, where the spacing is coarse, little 
difference is found in the temperature distribution. , For 
DELX = 0.001 cm, there is greater difference between the two 
because the equal spacing has limited the semi-infinite body to 
27 (0.001 cm) and fixes the temperature at the end point to 
970° C, thus not allowing the temperature to decline as rapidly as 
it should. For DELX - 0.0001 cm, the equal spacing method could 
not work, at all because 27 (0.0001 cm) is too small a range to 
define a semi -in finite body even for the small time constants under 
consideration. 

The second set of graphs, Fig. 3 (grid spacing DELX = 0.001), 
shows wide disparity between equal and unequal spacing whereas the 
first set of graphs (Fig. 2) (DELX = 0.01) shows good agreement. 
The smaller DELX needs more points to simulate the semi -infinite, 
one dimensional case and when restricted to IX = 28, fails to 
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970 r 


960 t“ 


950 - 


Initial belt Temp 


u 930 



900 f- 


0.002 0.004 

Position from Original Melt Surface* cm 

Fig. 3 Computed Temperature Distribution in the Semi-infinite Body at 
Different Times* t (initial grid spacing - 0.001 cm) 



Melt Temp 


970 



H 920 


910 i- 



Position from Original Melt Surface, cm 

Fig. 4 Computed Temperature Distribution in Semi-infinite Body at 
Different Times, t (initial grid spacing * 0.0001 cm) 



allow temperature away from the evaporating surface to decline 
rapidly because it is artificially pegged at x = 28 (0.001) to 
970°. The unequal spacing needs but six points to give equivalent 
extension and, when given 10 or 11 points, can adequately span a 
sufficient distance to simulate semi-infinity. At smaller 
DELX (0.0001) one cannot even attempt to use equal spacing with- 
out modifying the behavior at the last mesh point. For unequal 
spacing, 16 points will adequately represent the semi-infinite body 
for the times under consideration. 

Checking Program 

To check the program, GaAs single crystals will be grown 
with controlled dopant type, dopant concentration, growth rate, 
and temperature gradient, as shown in Table 5. The dopant con- 
centration, carrier mobilities, defect contents, ... will be mea- 
sured along several sections on each crystal. The results will be 
statistically analyzed and presented in future reports. 
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Crystal 

No. 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 
13 
16 

17 

18 

19 

20 


TABLE 5 GaAs CRYSTAL GROWTH SCHEDULE 

Temperature 


Dopant 

Concentr at ion 

Growth Rate 
in . /hr 

Gradient 

°C/in. 

Te 

1 

X 

10 17 

0.16 

8 

Si 

5 

X 

10 18 

0.22 

8 

Cr 

5 

X 

10 18 

0.28 

6 

Si 

5 

X 

10 17 

0.16 

6 

Zn 

5 

X 

io 18 

0.16 

4 

Cr 

1 

X 

10 18 

0.16 

10 

Te 

1 

X 

10 18 

0.22 

6 

Zn 

1 

X 

o 

oc 

0.28 

8 

Cr 

5 

X 

io 17 

0.10 

8 

Si 

1 

X 

10 18 

0.10 

4 

Si 

1 

X 

IO 17 

0.28 

10 

Zn 

1 

X 

io 17 

0.10 

6 

Te 

3 

X 

10 17 

0.28 

4 

Zn 

5 

X 

io 17 

0.22 

10 

Cr 

1 

X 

IO 17 

0.22 

4 

Te 

5 

X 

10 18 

0.10 

10 


To be grown after the results of above crystals 
are analyzed. 
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CONCLUSIONS 


We have developed a computer simulation program to study the 
phenomena of directional combined evaporation and solidification 
in binary alloys. A realistic phase diagram involving curved 
liquidus lines is used, and the program can be used for cases 
where the solid and liquid material parameters (e.g., specific 
heat, conductivity, diffusivity, ...) are functions of both tem- 
perature and solute concentration. The program works well if one 
follows the guidelines outlined in the report. The computed out- 
put results include the locations and velocities of movement of 
both the evaporation and solidification boundaries, and the tem- 
perature and concentration profiles in the semi-infinite alloy 
body at selected instants of time. 
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APPENDIX I 

IMPROVED COMPUTER PROGRAM 


^Kmnotoiim 


51 



FILE: BINCB6 FORTRAN Pi 


GRUMMAN 


DATA 


s y s i e t 


DIMENSION X (28 > ,T(28) ,TT (28) , C (28) ,CC{28) , 

*TEM(10) /A (84) 

NAMELIST /IN VAR/CEE, DS , DL ,TQ , CO , X K1 , R HO, GAMMA , RHOU , GAMMA U , EE, SI G 
*,ET,EA,EB,EC,ED,AU,BU, EMU , AV , EV,EMV,EK,IX 
D2 (X,F,V,G,Z,H) = ('(H-G)/(Z-Y)- (F-G)/(X-Y) )/{Z-X) *2. 

ABSl (X) =AMAX1 (1 . ,ABS(X) ) 

UE (V) = EK* (10.** (AU-BU/( 27 3, 12+V) ) } /SQRT (EMU* ( 2 73 . 12 ♦ V) ) 

V E ( V) = EK* { 1 A «** (AV-BV/(273.12+V) ) )/SQFT (EMV*(273. 12+V) ) 

FS (V) = (EB*(ET-V) +EA)*(ET-V) 

FL (V)= (ED* (ET-V) + EC}* (ET-V) 

XCL(V) =E T-2 , *V/ (EC+SQRT ((EC) **2 + 4.* (ED)* (V) ) ) 

XCSL (V) =ET+2.* V/ ( ( EC-EA) + SORT ( (E A- EG) ** 2+ 4 . * (EB-ED) * V) ) 

DFL (V) =- (2.* ED* (ET-V) + EC) 

11=1 

10=8 

R EAD (II, 100) IX,IAM,NIT,ItJ,NONCN 

100 FORMAT (1615) 

N ITH=N IT/2 
NITQ= NITH/2 
NITl=NITH+NITQ 

64 READ(II,101) ET,EA,EB,EC ,ED 
AQUAN=- (EA-EC) **2/ (4. * ( EB-ED) ) 

101 FORMAT (7E1C.0) 

READ (II, V 1) AU,BU, EMU, AV, EV,EM V,EK 

READ (II, I'M) CEE,DS,DL,TO ,C0, XKL,RH0, GAMMA , RHOU, GAMMAU,EE ,SIG 
ALS=XKL/(RHO*CEE) 

X KS= 1. 1* X KL 
A S S = 1 , 1* AL S/1 . 03 
A S= SQPT (ASS) 

AL=SQRT (ALS) 

READ (I I, 101) DELX,DELT,DEITM,DE1E,TF,S 
DEITK= AMAX1 (DELTM, DELT/1024.) 


BIN00010 
BIN00020 
BIN 000 30 
BIN00C40 
BIN00050 
BIN00060 
BIN00070 
BIN00C80 
BIN00090 
BIN00100 
BIN 00110 
BINOO 120 
BIN00130 
BIN 00 1 40 
BINOG 1 50 
BINOO 160 
BIN00170 

DTttnMan 

V AW W V 1 w V 

BIN 00 190 
BINO02G0 
BIN 00^10 
BIN00220 
BIN00230 
BIN 00240 
OIN00250 
BIN00260 
B1W00270 
BIN 00280 
BIN00290 
BIN003G0 
BIN 00 3 10 
BIN00320 


T2 = XCL (CG ) 

C0^ = FL (T 2) 

201 DO 1 1 = 1, IX 
IF (I- 2) 2,3,4 

2 X(1)=0. 

GO TO 1 

3 X (2) =DELX 
GO TO 1 

4 IF(IAM) 5,5,6 

5 X (I) =X (1-1) + DELX 
GO TO 1 

6 X (I) =X (1^ 1) +X (1-1) 

1 CONTINUE 

99 9 HR ITE (10, 107 )IX,IAM,NIT,IM, NONCN 
WRITE (10, 102) £T,EA , EB , EC, ED 
WRITE (10, 102) AU , BU , EMU , AV,BV,EMV,EK 
HR IT E( 10, 10 2) CEE,DS,DL, TO ,C0 ,XKL,RHQ, GAMMA 
*00° 

TPR=DELP 
R AT= 1 . 

T SI 1 = T0 
T 1 1 = TO 
CSL1=C0 


BIN00330 
BIN00340 
BIN00350 
BIN00360 
BIN0G370 
B IN 003 80 
BIN00390 
BING0400 
BIN00410 
BINQ0420 
BIN00430 
BIN00440 
BIN 00450 
BIN00460 
BIN 00470 
BIN00480 

,RHOU,GAMMAU,EE,SIG,T2,BIN00490 

BIN00500 
BIN00510 
BIN 00520 
BIN 00530 
BIN00540 
BIN00550 
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n. II M M A N 


DATA 


S Y S T £ n 



YI 1-0. 

ZI1-9. 

Ill =2 
112=2 
II 3 = 2 
D2T2=0. 

D2T3=0. 

D2TU=0. 

D2C3=0. 

D2C4=9. 

D2C2=0. 

D2T1 = '>. 

D2C1=^ . 

DTLDX=0. 

DTSDX=C . 

D2C5=H. 

D2T5=0. 

D2C 6=0 . 

D2T6=0. 

TIME =0. 

DCDXO =0 . 

DTD XO =0 . 

DELT S= DELT/8 • 

TIME1- TIJ1E + DELTS 
DO 1* 1= 1 , IM 

C (I) =C0 
CC { I ) =co 
TT ( I) =T0 

10 T (I) =T0 

I FL=0 

I FS=0 
111=112 

14 I T=0 

IF (IFS.EQ. 1) GO TO 1 99 
IF (TFL) 1 1 , 1 1, 20 

11 «0=OF(TO) 

yfl = VE (TO) 

IF 1= 1 

199 IF(IFS.NE.O) CSLl=CSl 
IF(IFS.EQ. 1) IFS=2 

DYDTO=UO*E«U*G SL1/RH0U+ VO*EMV* (1. -CSLl) /BHO 
HB0 = -EE*SIG*( 27 3. 1 2+TS II } **4-UQ* GAWflAU*CSLl -VO*GAMfl A* <1 . -CSLl) 
DCDTO = DCDXO* DYDTC- (UO-VO) *CSL 1 
DTDT / ' = DTDX n *D YDTO + HB^ 1 
20 YI2=YI1+DELTS* DYDTO 
IF <IFS.EQ .0) Z.I2 = YI2 
IF(IFS.NE.O) ZI2=ZI1+DELT5* DZDTO 
IF {ZI2 . GT. X (II 2*1) ) ZI2= (X (II 2) +YI2) /2. 

TSI2=TSI1+DELTS*DTDT0 
TI 2=TSI 2 

C SL2 = C SL 1 +DELTS*DCDTC 
IF(IFS.EQ.O) CL2=CSL2 
77 IF (IFS. EQ. 0) GO TO 777 
IF (IFL.GT.I)GO TO 877 

I I T=9 



BIN 00560 
BIN00570 
BIN 00560 
BINQ059U 
BIN 00600 
BIN 00 61 0 
BI NOG b 20 
BIN j0630 
BIN00640 
BIN 00650 
BIN0C660 
BIN 00670 
BIN00680 
BIN 00690 
BINU07G0 
BIN00710 
BINOO 7 20 
BINOO 73 0 
BIN OC 740 
BIN00750 
BIN 00760 
B1N00770 
B IN 0 0 7 6 0 
BIN00790 
BIN 00800 
BIN00810 
BI N 00620 
BIN 00830 
BIN00640 
BIN 00850 
UIN00860 
BIN 00870 
BIN00880 
BIN00890 
BI NOO 900 
BIN C091 0 
BIN 00920 
BIN00930 
BIN 0C940 
BIN00950 
BIN00960 
BIN00970 
BIN00980 
BIN00990 
BI N0 1 000 
B1N01010 
BINO 1020 
BIN01030 
BIN 01040 
BIN0 10 50 
BINQ 1060 
BIN01070 
BINO 1080 
BI N0 1 090 
BIN 0 1 100 
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PILE: BINCR6 FORTRAN Pi 


GRUMMAN 


DATA 


S X S T E t 


St 


877 CS2=FS(TI2) 

CL2 = PL(TI2) 

777 D2C2=D2 (ZI2,CL2,X (112) , CC (112) ,X (112*1) ,CC(II2+1)) 

CC (112)= (CC (II 2) +C ( 1X2) *»5*DELTS*(C2C1*D2C2) *DL > /2 . 

CALL MOTON (X (112) , CC(II2) ,ZI2,C12 ,X (112*1) ,CC(II2* 1) ) 

IF (ZI2.LT. X(II2-1) ) CC(II2-1) *CC(1I2) ♦ (X { II 2- 1) - X ( 1 12) ) * 

* (CL 2-CC ( II 2) )/ (ZI2-X(II2> ) 

IF (III -I 13-1) 83, 87, 84 

87 XP“ YI2 
CP=CSL2 
GO TO 184 

84 XP = X (I1 1-2) 

CP=CC (II 1 -2) 

184 D2C4 =D2 (XP ,CP, X (III- 1) ,CC (I1 1-1) ,ZI2,CS2) 

CC(III-I) ®C(II1-1) +.5*DELIS»(D2C3+C2C4)*DS 
CALL MOTON (X(II1-1) ,CC(II1-1) ,XP,CP,ZI2,CS2) 

83 IF (II2.EQ.II1.0R.ZI2.LT.X(II2-1) ) GC TO 85 
IF (II2-II3.GT. 1) GO TO 185 
XP=YI2 
CP=CSL2 
GO TO 51 

185 XP=X (II2-2) 

CP=CC (II2-2) 

51 CC (II2-1)=CP*(X(II2-1) -XP)*(CS2 -CP)/ 

*(ZT2 -XP) 

8 5 D2T2=D2 (ZI2,TI2,X (II 2) ,TT (II 2) ,X (II 2* 1) ,TT (112*1) ) 

IF (D2T2.GT.0.) D2T2=0. 

TT (II 2) - (TT (112) *T (112) *. 5*DELTS* (E2T 1*D2T 2) *ALS )/2. 
CALL MOTON ( X (1 12) ,TT (112) ,ZI2,TI2,X(II2+1) ,TT(II2*1)) 

IF (IF S • EQ • 0) GO TO 485 

48 5 IF (ZI2.LT. X{ II 2-1) ) TT (1 12 -1 ) =TT ( 1 1 2) + (X(II2-1) -X(II2)) * 

* (TI2-TT (112) ) / (ZI2-X (II 2) ) 

IF ( II 1-II3- 1) 69,1 69,269 

169 TP =TSI 2 
XP= YI2 
GO TO 16 

269 TP=TT (II 1 -2) 

XP=X (I1 1- 2) 

16 D2TU = D2 (XP,TP, X (I1 1-1) , TT (III -1 ) ,ZI2,TI2) 

TT (III — 1 ) - T (1 11-1) *.5*DELTS*(E2T3*E2T4) *ASS 
CALL MOTON (X (I1 1-1) ,TT ( 1 1 1 - 1 ) ,XP,TP ,ZI2,TI2) 

69 IF(II2.EQ.II1.0R.ZI2.LT.X(II2-1) ) GC TO 86 

52 IF (II2.LT. II3-1) GO TO 186 
XP= YI 2 

T P=TSI 2 
GO TO 352 

186 XP =X ( 112- 2) 

TP- TT ( II2-2) 

352 TT (II 2- 1 ) = TP * (X (I I 2- 1 ) -X P) * ( TI 2 -TP)/ 

* ( ZI2 -XP) 

86 IF (IPS. EQ.O) GO TO 299 

DCLDX= (CL2-CC (112) ) / (ZI2 - X (II 2) ) 

IF (D2C2.GT.0.) DCLDX=DCLEX-D2C2* (X (112) -ZI2) /2 . 

IF (D2C2.LT. 0.) D2C2=9, 

DTLD X= (TI2 -TT ( 112) ) / (Z I 2-X ( 112} ) 


SI N0 1 1 10 
BIN01 120 
BIN01130 
BIN 01 140 
BING 11 50 
BINO1160 
BIN01170 
BIN 01160 
BIN01190 
BINO 1200 
BIN01210 
BIN01220 
BINO 1230 
BIN01240 
BIN01 250 
BINO 1260 
BIN01270 
BINO 1280 
BINO 1290 
BIN01 300 
BIN 01310 
BIN 0 13 20 
BINOl 330 
BINO 1340 
BIN01350 
BIN01 360 
BIN 01 370 
BINO 1 380 
BIN01 390 
BIH01400 
BIN01410 
BINOl 420 
BIN 01 430 
BINO 1440 
BINO 1450 
BIN 01460 
BINO 1470 
BINO 1 480 
B1N01490 
BIN0 1500 
BINOl 510 
BINO 1520 
BING 1530 
BINO 1540 
BIN01550 
BIN 01 560 
BINOl 570 
BINOl 580 
BI N0 1 590 
BINOl 600 
BIN01610 
BINOl 620 
BINOl 630 
BIN01640 
BINOl 650 
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PILE: BINCR6 


FORTRAN PI 


GRUMMAN 


DATA 


SYSTEM 


IF (D2T2 .LT.O.) DTLD X=DTLD X- , 5*D2T2* {X (112) -ZI2) 

IF (II3.EQ.II1) GO TO 386 

XP=X(II1-1) 

TP = TT ( I1 1- 1) 

CP=CC (II1-1) 

GO TO 486 
386 X P =YI2 
TP=T5I 2 
CP=CSL2 

DC SDX= (CS2-CP) /(ZI2-XP) 

DTSD X= (TI2-TP) /(ZI2-XP) 

GO TO 686 

486 DCSDX= (CP-CS2)/(XP-ZI2) -D2C4* (XP-ZI 2) /2. 

DTSD X= (TP-TI2) /(XP-ZI2) -D2T4* (XP-ZI 2) / 2. 

686 DZDT= (DL*DCLDX-DS*DCSDX) / (CS2-CL2) 

DZDTT= ( X K S * DT S DX-X K L* DT LDX) / (RHC* GAMMA) 
FSL=DZDT* (CS2-CL2) /DZDTT 
IF (FSL.GT.C. .OS.FSL.LT. AQUAN ) GO TO 772 


6IN01660 
BIN 016 70 
BI NO 1660 
BIN01690 
bINO 1700 
BIN01710 
BIN0 1720 
31NO1730 
BINO 1740 
UIN0 17 50 
BINO 1 760 
BIN01770 
BIN01780 
BING1790 
BIN 0 1 800 
BINO 1810 
BIN01820 
B IN 0 1 8 30 


TII= XCSL (FS L) 

GO TO 771 

772 TII-ET 

771 COE1=XKS/ (ZI2-XP) 

COE2 = XKL/( X ( II 2) - ZI2) 

5 86 TI = (RHQ*GAMMA*DZDT+COEl* (TP-D2T4*. 5* (XP-ZI 2) ** 2 ) +COE2* (TT(II2) - .5 
**D2T2* (X (112)- ZI2) ** 2 ) ) / (COEl +COE2) 

773 IF (TI.LT.TSI2 .AND. TI.LT. Til) TI=TII 

IF (TI .LT.TSI2 .OR. TI ,GT.TT( II 5+1) ) TI=TT ( I 12) 

IF (DZDT. LT.O. .AND. DZDTT, GT.O.) DZDT=DZDTT 
IF (ABS (TI-TI2) -1. £-5*ABSl (TI+TI2) ) 5 6 7,587,770 
587 IF (ABS (DZDT- DZDTl ) - 1.E- 3* ABS1(DZDT*DZ DTI)) 2 98,298,770 
770 TI 2 ~ (TI+TI 2) /2 . 

IF (DZDT .LT.O.) DZDT = DZDT 1 
DZDT1= (DZDT ♦ DZDTl) /2. 

ZI 2 =Z1 1 + « 5 *DELTS* ( DZDTl ♦DZDTO) 

IF (7.12 . GT. X (II 2 + 1) ) ZI2= (X (112) +YI2) /2. 

IIT-IIT+1 

IF (IIT.GT. 5) GO TO 298 
GO TO 877 * 

298 IF(TI.GT.TT(II2) ) T I=TT (II 2) 

IF (TI.LT. TSI2) TI=TS 12 

IF ( II 1-II 3-1) 398,498,598 

598 D2T6 = D2 (YI2,TSI2,X (113) ,TT(II3) , X ( II 3+ 1 ) , TT <113+ 1) ) 

D2C6-D2 (YI2,CSL2,X (II 3) ,CC(II3) ,X(II3+1) ,CC(II3+1)) 

DTDX1= (TSI2-TT (1 1 3) ) / ( YI2 -X (1 13) ) -D2T6* (X(II3) -YI2) /2. 

DCDX1= (CSL2-CC (1 13) )/ (Y I2-X ( 1 13) ) -D2C6* (X ( 1 13) -YI2) /2. 

GO TO 599 

498 DTDX1 = DTSDX + D2T4* ( X P-Z I 2) 

DC BX1 = DCSD X + D2C4* ( X P-ZI 2) 

D2T 6=D 2T4 
D2C6=D2C4 
GO TO 599 
398 DCDX1=DCSDX 
DTDXI^DTSDX 
GO TO 599 

299 DTDXl - (TT (II 2) -TI2)/(X (II 2) -YI2) 5*D2T2* (X (II 2) -II 2) 


BIN01640 

BIN01850 

BIN01860 

BIN01870 

BIN 01 880 

BIN01890 

BIN01 900 

BIN 019 10 

BI N0 1 920 

BINO 1930 

BIN0 19 40 

BINO 1 950 

BIN01960 

BIN01 970 

BIN 01 980 

BINO 1990 

BIN02000 

BIN 020 10 

BIN020 20 

BIN02030 

BIN02040 

BIN02050 

BIN02060 

BIN02070 

BIN02080 

BIN02090 

BIN02100 

BIN 021 10 

B1N02120 

BIN02130 

BIN 02 140 

BIN02150 

BIN02160 

BINC2170 

BIN 02 1 80 

BIN02190 

BIN02200 


ORIGINAL PAGE IS 
OE POOR QUALM 
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PILE: BINCR6 


FORTRAN PI 


G R U M tl A N 


DATA 


SISTER 


DCDX1= {CC (112) -CL2 ) / (X ( 112) -YI2) - . 5*D2C2* ( X (112) -YI2) 

59 9 01=0E(TSI2) 

VI = VE {TSI 2) 

HB 1— ** EE*S IG* (273.12 + TSI2) **4-Ul*GAMHAU*CSL2-V1 *GAMttA* (1 .-CSL2) 

DYDT1=U1*EMU*CSL2/RHOU+V1*ISV* ( 1 . -CSL2) /RHO 

DCDTl =DC D X 1* DY DTI" (Ul-Vl) *CSL2 

DTDTl-DTD X 1* DYDTl + H Bl 

IF (IFL.LT.2) GO TO 76 

IF ( IF S . GE. 1) GO TO 174 

IF ( TSI2 . GT. TT ( 112) ) GO TO 399 

TT2= XC L (CSL2) 

IP(TSI2.GT.TT2) GO TO 174 
IF (TSI2.GE.TT2". 05) GO TO 400 
IF (DELTS. LE.DELTK) GO TO 400 
GO TO 399 

76 YT=YI 1+ . 5* (DYDT9 + DYDT1) *DELTS 
IF(IPS.EQ.O) ZI = YI 
TSI =TSI1+.5*( DTDTO + DTDT1) *DEITS 
IF (TSI. LT.O .) TSI=TSI 1 
CSl=CSLl+. 5* (DCDT0+ DCDTl) * DELIS 
IF (IFS. EQ.O) GO TO 73 
IF (IIT.GT » 5) GO TO 7 0 

IF (TSI.GT. TT (112) . OR. TSI. GT.TI2) GO TO 70 
7 3 IF (ABS (YI-YI2) -1.E-6*ABS1{YI+YI2)) 74,74,70 

74 IF (ABS (CSL-CSL2) -1 . E-4*AB Si (C SL+CSL2) ) 75,75,70 

75 IF (ABS (TSI "TSI 2) -1 , E-5* ABS1(TSI+TSI2))7,7,7C 
70 TSI2= (TSI + TSI2) /2. 

IF (IFS. Eg. 0) TI 2=TS 12 
YI2= { YI+YI2) /2 . 

IF (IFS.EQ.O) ZI2=YI2 
C SI2= (CSL+CS L2) /2. 

IF (IFS.EQ.O) CL2=CSL2 

IF(ZI2.LT. (1 .-S) *X(II2) +S*X(II2-1) > GO TO 24 
22 17(111-112) 24,96,96 
9 6 112=112+1 

4 6 D2C1=D2(ZI1,CL1,X (1 12) ,C(II2) ,X (1 12+1) ,C(II2 + 1) ) 

D2T1 = D 2 (ZI1,TI1,X (I I 2) ,T(II2) ,X(II2+1) ,T{II2+1) ) 

24 IF (IT-NITQ) 48,174,160 
16 ^ IF (IT-NITH) 48,47,48 

47 IF (DEL TS-DELTK ) 48,48, 53 

53 DELTS=DELTS/RAT 

152 IF (RAT-1. ) 1 53,1 53, 154 

153 RAT=2.*RAT 

DELTS= DELTS/2 . 

GO TO 152 

154 TI«E1=TIME+DELTS 
212 DO 45 I=II 3, IH 

TT (I) = (TT ( I) +(RAT-1.)*T (I)) /RAT 
4 5 CC (I)= (CC (I) ♦ ( RAT- 1.)*C(I)) /RAT 
148 R A T = 1 , 

48 I T=IT+ 1 

IF ( IT. EQ. N ITH+ 1) GO TO 20 
IF (IT-NITL) 161,174,161 
161 IF (IT-NIT) 77,77,402 
402 IF(DELTS.GT.DELTK) GO TO 399 
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BIN02500 
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BIN02620 

BIN 02630 
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BIN02650 
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PILE 


BINCR6 FORTRAN Pi 


GRUMMAN 


DATA 


SYSTEM 


IFI=2 
GO TO 77 

26 IF ( NONC N ) 99,48,40 
7 IFI=IFL+1 
DZDT1 =DZDT 

IF(IFS.NE.O) ZI2=ZI1+.5*DELTS* (DZDT0+DZDT1 ) 

TSI2=TSI 
YT2 =YI 

IF (IFS.EQ.9) ZI 2 = YI2 
CSL2=CSL 

IF (IFS.NE.O) TI2=TI 
IF (TI2.LT.TSI2) TI2=TSI2 
IF (IFS.EQ.O) TI2=TSI2 
IF(IFS.EQ.O) CL 2= CS L2 
GO TO 77 

399 I F (DELTS . LE. DELTK) GO TO 26 
IT= NITH 
I FL = 1 
GO TO 153 

117 IF (TIME. £Q. TIME1) TI M21 = TI ME1 +DELTS 

TINE=TIME1 

IF (R AT. NE.2. . AND. RAT. NE.O.) CELT S=DELTS/R AT 
R AT=1 . 

IF(IT-NITQ) 82,82,81 
82 IF (DELTS. GT. DELT/2 . ) GO TC 81 
DELTS=DELTS+DELTS 
R AT=2 . 

81 TI ME 1=TIME1+ DELTS 
282 D2C1-D2C2 
D2C3=D2C4 
D2C5=D2C6 
D2T3=D2T4 
D2T1=D2T2 
D2T5=D2T6 
D YDTO = DYDT 1 
DCDT9=DCDT 1 

uo=ui 

V0=V1 

DTSDX1=DTSDX 

DCSDX1=DCSDX 

IF(IFS.NE.O) DC DXO = DCDX1 

IF (IFS.NE.O) DTDI9= DTDX 1 

IF(IFS.NE.O) DZDT0 = DZDT1 

CSLl=CSL2 

TSI1-TSI2 

YI 1=YI 2 

ZI 1 = ZI 2 

IF (IFS.NE.O) CS1=CS2 

HB0=HB1 

CL1=CL2 

TI1=TI2 

IF (YI2.LT. X(II3) ) GO TO 410 
113=113+1 

D2T5=D2(YI1,TSI1 ,X(II3) ,TT(II3) ,X(II3+1) ,TT(II3+1)) 
D2C 5= D2 ( Y 1 1 / CSLT, X (113) ,CC (II 3) ,X(II3+1) ,CC (II3 + 1) ) 
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VILE: BINCR6 FORTRAN Pi 

410 IF (III .EQ. 112. OR. ZI2.LT. S* 
111=112 

IF (II1-II3-1) 33,310 , 110 
310 XP = Y1 1 
CP=CSL1 
TP=T5I 1 
GO TO 210 
110 xP=X(II1-2) 

CP=CC(II1-2) 

TP-TT (III -2) 

210 D2C3=f)2 (XP,CP, X{111-1) ,CC < 
D2T3=D2(XP,TP,X(II1-1> ,TT< 
IF (II1-II3.ME. 1) GO TO 33 
D 2T 5=D 2T3 
D2C5=D2C3 
GO TO 33 
174 111=112 

IF 'ZI 2 . LT. X {II 2-1) ) III=Ii: 
IF (II1-II3.LT. 2) 00 TO 1 ! 

TF (II1-II3.GT. 2} GO TO 29 
TT (113 )= (T (113) /DELTS+ . 5*J 
*+TT (113 + 1 ) *ASS/(X(II3+1)-; 
* (113) -YI2 )/(X(II3+1) -Y I2|] 
*+1)-X (II 3) ) + 1./(X(II3> -YI i 
CC (TI3) ={C (113) /DEL TS + . 5*1 
*+CC (II3+1)*DS/ (X (II 3 + 1) -X 
*1 3) - YI2) / (X(II 3 + 1) -YI2) )/ 
*-X (113) ) +1./{X (TI3) -YI2)) ; 
GO TO 18 

29 X{II3-1) = (YI2+YI1) /2 * 

C (TI3-1) =CSL 1 
CC(II3-1) =CSL2 

T (IT3-1) -TS1 1 
TT (II3-1 } =TS 12 
CALL TRIST(X,T,TT,II3-1 ,1 
CALL TRIST ( X, C ,CC, 113- 1, 
GO TO 18 

27 IF(T(IH-1) -T0+1.E0) 30, 

30 IF (IM-IX) 15,32,32 
15 IM=If1+1 

T ( I M) = TC 
C(IM) =CC 

32 TT(IM) =T(IM) 

CC (IM) = C<IM) 

GO TO 14 

18 IF (DELTS.LT. DELTK) DELTS=D 
CALL TRIST(X,T,TT,III, Ifl, 
CALL TRIST (X,C,CC, III, IM, 
DO 21 1= 112 ,IM 
IF (I.EQ. II 2) GO TO 21 
21 CONTINUE 

219 IF (IFL.EQ.2) GO TO 117 
GO TO 48 

33 IF (TIME -TPR) 50,34,34 
50 IF (NONCN) 90,98,54 


GRUMMAN DATA SYSTEM 


*X ( 1 1 2} + ( 1.-S) *X(II2-1) ) GO 10 33 


(II1-1) ,ZI 1 , CL 1) 
(III -1) ,ZI1,TI1) 


A SS* D 2T 5 

X (II 3) }/ (X (II 3 + 1) -YI2) +TSI2*ASS/ 
)/(1./DELTS+ASS/(X (II 3+ 1) - YI2) * ( 
2 ))) 

E5*D2C5 

(113) ) / (X ( II3+1) - YI2) +C SL2*D $/ (X 
(1. /D ELI S + DS/(X (1X3+1) -YI2) *(1./ 
) 


II- 1 , ASS , DELT S , A) 

III- 1, DS, CELTS, A) 

32,32 


ELTK ■ 

ALS , DELIS , A) 
DL, DELTS, A) 


31 N 033 1 0 
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FILE: BINCR6 FORTRAN Pi 


G B U H H A N 


DATA 


SYSTEM 


54 

42 


143 

^ f\ o 

98 


IF(IT-NIT) 98,98,42 
IF (DELP.EQ.G.) GO TO 42 
TPR=TPR+DELP 

WR ITE ( 10 . 100) IT f IM ,111,112 ,113 

WRITE (10, 102) TIME, Y II ,C SLl , TS 1 1 , Ul , V 1 ,H Bl , DYDTl , DT DT 1 , DC DT 1 , DELTS 
IF (IFS . NE .0) WRIT E ( 10 , 1 0 2) Z1 1 , CS 1 ,CL 1 ,T1 1 , DZDTl 
WR ITE (10, 10 2) (TT (I) ,I = H 3 , 1 M } , (CC (I) ,I = II3,Ifl) 

FORMAT ( 5E1 4 . 6) 

IF (TIME -TF) 98,99,99 

DO 97 1=113 , IM 

TTT=TT (I) + RAT* (TT (I) -T ( I) ) 


T (I) =TT (I) 
i Ti(i) =TTT 
97 CONTINUE 

DO 197 I=II3,IH 
CCC=CC (I) ■♦■BAT* (CC (I) -C (I) ) 

C ( I) =CC (I) 

197 CC(I)=CCC 
IFL-1 
GC TO 27 
4*^ CSS=FS (TI2) 

CLL - FL (T 1 2) ^ 

DEL7= (CLL-CL2) / ( (CLL+CC (112) ) - 2. *CSS) * ( X (II 2) - YI 2) 

ZI2=YI2+DELZ 
DZDTf =DELZ/DELTS 
DZDT1=DZDT0 
CL2=CLL 


C S2 = CSS 
DCDXO =0 . 

PTDXO =0 • 

L * *" D X 1 = 0 * 

WPITE(IO,103)'TIME1,YI2,TI2,ZI2,CS2,CL2,DELZ,DELTS,DZDTO 

^3 FOPMAIC SOLIDIFICATION HAS BEGUN */(5El4.6)) 

I F S -1 
GO TO 174 
99 RFAD(II, INVAR) 

IF (IX. GT. 0) GO TO 999 

STOP 

END 
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APPENDIX II 


IX 

IAM 

NIT 

IM 

NONCN 


IT 

NITH 

NITQ 

NITL 

IFL 

IFS 


111 

112 

113 

1 14 
III 
DELT 
DELTM 
DELTK 
DELTS 
TIME 
TIME1 
DELP 
TPR 
TF 


GLOSSARY OF PROGRAM PARAMETERS 


- maximum number of points in mesh, <28 

- spacing option: 0 indicates equal, 1 unequal doubling 

= maximum iteration count 

= number of points in mesh 

= nonconvergence option: 1 indicates proceed and printout 

.0 indicates proceed but do not 
printout 

-1 indicates program stop 

= actual iteration count 
= half of NIT 
= quarter of NIT 
= 3/4 NIT 

= indicator of convergence: 2 on convergence, 

< 2 before convergence 

= indicator of beginning of solidification: IFS - 0 before 

solidification, 


IFS > 0 after 
solidification 

= grid point reference for solid side of mesh 
= grid point reference for liquid side of mesh 
= grid point reference for evaporation boundary 
- grid point reference for point after evaporation boundary 
= grid point reference for point after solidification boundary 
*» maximum time interval (step size) 

= minimum time interval 

= larger of quantities DELT/ 1024 and DELTM 
= current time interval 
= time at beginning of time interval 

= time at end of time interval _ BLANK H0T FILM® 

= v 'time print interval 


fBECEDING PAGE 


= time for printing results 
= final time 
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YI, YI1, YI2 
ZI1, 212 

TI, Til, Til, TI2 
TSI, TSI1, TSI2 
CSL, CSL1, CSL2 
CS1, CS2 

CL1, CL2 

DZDT, DZDTT, DZDTO, 
DYDTO, DYOT1 
DTDTO, DTDT1 
DCDTO, DC DTI 

DTDXO, DTDXl 

DTSDX 

DTLDX 

DCDXO, DCDXl 
DCSDX 

DCLDX 
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= values of y (evaporation boundary) 

= values of z (solid-liquid boundary) 

“ temperatures at solid-liquid boundary 

“ temperatures at evaporation boundary 

= concentration at evaporation boundary 

= concentration of solid at solid-liquid 
boundary 

= concentration of liquid at solid-liquid 
boundary 

dz 

DZDT1 = — derivative of solid-liquid boundary 

dy 

= derivative of evaporation boundary 


dT 

= derivative of temperature at evapo- 
ration boundary 

= ^“r derivative of concentration at 
evaporation boundary 



partial derivative for tempera- 
ture at evaporation boundary 


/dT s 

-(sr) 


partial derivative of tempera- 
ture in solid at boundary 



partial derivative of tempera- 
ture in liquid at boundary 

partial derivative of concentra- 
tion at evaporation boundary 


-(&■): 


partial derivative of concentra- 
tion in solid at solid-liquid 
boundary 


y z partial derivative of concentra- 
tion in liquid at solid-liquid 
boundary 



D2T1, D2T2 


D2T3, D2T4 
D2T5, DdT6 

D2C1, D2C2 

D2C3, D2C4 

D2C3, D2C4 

D2C5, D2C6 

UO, U1 
VO, VI 
HBO, HBl 

EMU, EMV 

RHOU, RHO 
GAMMAU, GAMMA 
AU, BU, AV, BV, EK 

UE, VE 

ET, EA, EB, EC, ED 
FS, FL • 


^ 2 _ 
,d T, 


(rr) 

ox 


*2_ 
, o T 


(rr) 

ox 


(ft 

dx y 


a 2 c. 


(-/) 

dx 2 


* 2 „ 

,o C 


(pr) 

ox z 


b 2 c 


(!% 

ox y 


at liquid side of solid-liquid 
boundary 


at solid side of solid-liquid 
boundary 

at evaporation boundary 


at liquid side of solid-liquid 
boundary 


at liquid side of solid-liquid 
boundary 


at solid side of solid-liquid 
boundary 

at evaporation boundary 


- rates of evaporation of solute 

= rates of evaporation of solvent 

“ heat balance sum of evaporation and 
radiation terms 


= molecular weight of solute and solvent 
atoms 

= density of solute and solvent 

- specific heats of solute and solvent 

- evaporation constants for solute and 
solvent 


= arithmetic function definition for 
evaporation rates 

= phase diagram constants 

= arithmetic functions for solidus and 
liquidus curves 
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AS, AL 


temperature diffusion coefficients 

ASS, ALS 


squares of temperature diffusion 
coefficients 

XKS, XKL 

= 

k > for interphase boundary equa- 

tion 

DS, DL 


mass (concentration) diffusion coeffi 
cient 

EE, SI G 


radiation constants e, a 

TO, CO 


initial temperature and concentration 
distribution 

COO 

= 

equals CO 
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APPENDIX III 


LIST OF SYMBOLS 


A ,B 
u u 

V B v 

c 

C 

C i 


e a' e b 

e c ,e d 


K 


M 


u 


M 


v 


k,/p c 
£ v 

k /p c 
S V 

evaporating constants for solute 

evaporating constants for solvent 

specific heat, cal/g/°C 

solute concentration 

solute concentration in liquid 

initial solute concentration in liquid 

solute concentration in solid 

3 

density of solid, g/cm 

3 

density of liquid, g/cm 
melting point of solvent, °C 
coefficients of solidus 
coefficients of liquidus 
liquidus equation 
solidus equation 

constant in evaporation equations - 5.83 x 10 
thermal conductivity of liquid, cal/sq cm/ cm/sec/ °C 
thermal conductivity of solid, cal/sq cm/cm/sec/°C 

molecular weight of solute 
molecular weight of solvent 


-5 



time constant 
time, second 
temperature, °C 

temperature at solid-liquid interface, °C 
temperature in liquid, °C 
initial melt temperature, °C 

temperature in solid, or surface temperature, °C 
surface temperature, °C 
final surface temperature, °C 

2 

solute evaporating rate, raol/cm /sec 

2 

solvent evaporating rate, mol/cm /sec 

position and temperature -dependent variable 

distance from initial melt surface, cm 

distance at phase change boundary, cm 

rate of movement of phase change boundary, cm/ sec 

distance at phase change boundary, cm 

constant 

emissivity coefficient 
3 

density, g/cm 

3 

density of solute, g/cm 

3 

density of solvent, g/cm 
latent heat of fusion, cal/g 

latent heat of fusion, or specific heat; of solute, 
cal/g or cal/g/ °C 

latent heat of fusion, or specific heat, of solvent 
cal/g or cal/g/ °C 

Stefan-Boltzmann constant = 1.35 x 10 
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